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Twisted complex superfluids in optical lattices 
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We show that correlated pair tunneling drives a phase transition to a twisted superfluid with a complex order parameter. 
This unconventional superfluid phase spontaneously breaks the time-reversal symmetry and is characterized by a 
twisting of the complex phase angle between adjacent lattice sites. We discuss the entire phase diagram of the extended 
Bose-Hubbard model for a honeycomb optical lattice showing a multitude of quantum phases including twisted 
superfluids, pair superfluids, supersolids and twisted supersolids. Furthermore, we show that the nearest-neighbor 
interactions lead to a spontaneous breaking of the inversion symmetry of the lattice and give rise to dimerized density- 
wave insulators, where particles are delocalized on dimers. For two components, we And twisted superfluid phases with 
strong correlations between the species already for surprisingly small pair-tunneling amplitudes. Interestingly, this 
ground state shows an inflnite degeneracy ranging continuously from a supersolid to a twisted superfluid. 


The time-reversal symmetry of the Schrodinger equation al¬ 
lows us to describe superfluid ground states with real wave 
functions in accordance with Feynman’s no-node theorem^ 
This principle holds unless the Hamiltonian breaks this sym¬ 
metry either explicitly or spontaneously, which is typically 
the case for higher orbitals or spin-orbit coupling^. Re¬ 
cently, experiments with binary mixtures of ultracold bosonic 
atoms in a honeycomb optical lattice have observed a multi¬ 
orbital superfluid^. Most strikingly, time-of-flight measure¬ 
ments have revealed a complex superfluid order parameter. 
The complex phase angle of the superfluid order parameter 
twists between neighboring lattice sites without accumulating 
a total flux. Furthermore, a phase transition between a nor¬ 
mal and this so-called twisted superfluid was observed. The 
origin of this unconventional superfluid phase lacks a con¬ 
clusive theoretical understanding. In fact, previous theoreti¬ 
cal studies"^’^ have found no indication for a transition to the 
twisted superfluid phase. Here, we show that correlated pair 
tunneling can drive the quantum phase transition to the twisted 
superfluid ground state. This quantum phase is remarkable in 
several ways. It spontaneously breaks the time-reversal sym¬ 
metry, combines off-diagonal long-range order with diagonal 
short-range order and is crucially based on beyond-mean-fleld 
correlations. Interestingly, the twisted multi-orbital superfluid 
has a strong connection to mechanisms discussed for super¬ 
conducting materials. Two-orbital order parameters have also 
been experimentally observed^"^ in high-temperature super¬ 
conductors and play a central role for the understanding of this 
phenomenon. A spontaneous time-reversal symmetry break¬ 
ing associated with complex order parameters has been ob¬ 
served in the pseudo-gap phase of the high-temperature su¬ 
perconductor Bi-2212 with ARPES^, as well as in Sr 2 Ru 04 
using a muon spin relaxation measurement^^ and UPts via the 
polar Kerr effect^In this respect, the twisted superfluid¬ 
ity interlinks many-body quantum gas systems with supercon¬ 
ducting materials. 

The Hubbard model is a primary description for strongly 
correlated electrons in lattices. It accounts for the single¬ 
particle tunneling between neighboring sites with amplitude 
J and the interaction t/ of a pair of particles on the same lat¬ 
tice site. Its counterpart for bosonic particles is the Bose- 


Hubbard model. For weak interactions the ground state 
is a superfluid phase (SF), whereas for strong interactions 
a localization of the particles is favored forming the Mott 
insulator phase (MI)^"^"^^. Theoretically, extended Bose- 
Hubbard models with off-site interactions^^ have been studied 
predicting pair superfluid phases charge-density-wave 

insulatorsand supersolid phasesThe lat¬ 
ter two phases exhibit a density modulation with respect to 
neighboring lattice sites and break the inversion symmetry of 
the lattice Hamiltonian spontaneously in the thermodynamic 

limit26,27 

We study an extended Bose-Hubbard model with nearest- 
neighbor interaction and pair tunneling on a honeycomb 
lattice using cluster-Gutzwiller theory. This method takes 
into account short-range correlations exactly and simulta¬ 
neously describes the thermodynamic limit of inflnite two- 
dimensional lattices. The extended Hubbard processes give 
rise to several novel correlated quantum phases. First, we 
discuss the driving mechanism behind the twisted superfluid 
phase and the spontaneous breaking of the time-reversal sym¬ 
metry. Second, we show how nearest-neighbor interaction 
supports dimerized density waves, where the particles are de¬ 
localized on dimers and the translational symmetry is broken. 
This insulator is surrounded by a dimerized supersolid ex¬ 
hibiting the same density modulation. Third, we extend our 
study to a quantum gas mixture of two spin components. Re¬ 
markably, the two-species twisted superfluidity persists even 
at small pair-tunneling amplitudes allowing its experimental 
realization with ultracold atoms^. The twisted superfluid is 
characterized by strong correlations between the two compo¬ 
nents and exhibits an inflnite ground-state degeneracy. 

Results 

We study the extended Bose-Hubbard model with two addi¬ 
tional processes arising from the interaction of particles on 
adjacent sites known as nearest-neighboring interaction V and 
(correlated) pair tunneling P (see Fig. la). As in the experi¬ 
mental realization^, we consider a honeycomb lattice geome¬ 
try. First, we restrict our calculations to bosonic particles with 
a frozen spin degree of freedom. For a single component, the 
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FIG. 1: Phase diagram of the single-component extended Bose-Hubbard model, a Schematics for the processes in the extended Bose- 
Hubbard model, b Critical value of J/U for the transition to the superfluid (or supersolid) phase as a function of the effective chemical potential 
lieslU and the off-site interaction amplitudes P/U = V/U of pair tunneling and nearest-neighbor interaction. The three-dimensional phase 
diagram is obtained using a two-site cluster. The colors indicate the distinct quantum phases. Insulator phases with integer site occupancy are 
the Mott insulator (MI, dark green), c the charge-density-wave insulator (CDW, light green), and d the maximally-imbalanced CDW insulator 
(CDW2, gray). The phases e-h have a non-vanishing superfluid order parameter (indicated by a blur), e The twisted superfluid (TSF, yellow) 
is characterized by a complex order parameter (d^), where the arrows indicate the complex phase angle arg(di). f The twisted supersolid 
(TSS, orange) phase has an additional density wave, g Supersolids (SS) combine a real order parameter with a density wave (not shown in b). 
b Dimerized supersolids (DSS) are characterized by a long-range density wave with no density gradient within unit cells, i At quarter integer 
Ailing, dimerized density-wave insulators (DDW, blue) appear where particles are delocalized on two sites, j In addition, dimerized insulators 
with other fractional flllings such as 5/18 are possible. For large values of P and /x, pair superfluids (PSF, red) and pair supersolids (PSS, 
brown) can be found. 


Hamiltonian of the extended Bose-Hubbard model reads 
H =^^hi{hi - 1) - J 

i (ij) i 

p ( 1 ) 

+y hihj + 

where the operators hi (dt) annihilate (create) a particle on 
site i and hi = dtd^. The brackets (x, j) indicate pairs of 
nearest neighbors i and j. The first line of the Hamiltonian 
is the standard Bose-Hubbard model with on-site interaction 
energy U, tunneling amplitude J and chemical potential /x. 
The second line describes off-site interaction processes V and 
P giving rise to the exceedingly rich quantum phase diagram 
shown in Fig. lb. For neutral ultracold atoms, it is valid to 
assume a contact interaction potential. As a consequence, the 
amplitudes of the nearest-neighbor interactions V and the pair 
tunneling P are equal. 

The phase-diagram in Fig. lb is obtained by means of 
a cluster-mean field theory. Within each cluster quantum- 
mechanical correlations are treated exactly whereas the 
boundary is coupled self-consistently to the mean field. Un¬ 
like in the standard Hubbard model, the respective factoriza¬ 
tion of the Hamiltonian (1) leads to three distinct mean fields, 


namely (hi), {hi) and (hihi) (see Methods). The increased 
complexity limits the numerical treatment to clusters of up 
to four sites. The phase diagram Fig. lb covers the whole 
three-dimensional parameter space of the Hamiltonian (1) for 
P = V. For small ratios J/U insulator phases dominate the 
phase diagram except for the pair superfiuid (PSF)^^“^^ and 
pair supersolid where the order parameter {hi) 

vanishes but the pair order parameter {hihi) takes a finite 
value. 

Twisted complex superfluidity 

The most striking feature is the appearance of the twisted su¬ 
perfluid phase (TSF). As illustrated in Fig. le, this unconven¬ 
tional superfiuid is characterized by a twist of the complex 
phase angle of the superfiuid order parameter between neigh¬ 
boring sites, whereas the density {hi) is homogeneous. Ex¬ 
pressing the local order parameter as {hi) = (j)iU^\ the differ¬ 
ence of the complex phase angle of adjacent sites is ^ = Oi—Oj . 
The ground state breaks the time-reversal symmetry and is 
two-fold degenerate with a twisting of either 0 or —0. The 
two degenerate many-particle states {t/, T^p} are connected 
by the time-reversal symmetry operation Tpj = t/*. In Ref. 3, 
the complex superfiuid ground state has been observed ex¬ 
perimentally. It has been linked to the emergence of a two- 
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mode superfluid order parameter on the basis of a simple vari¬ 
ational mean-fleld model. Previous theoretical studies have 
found only real-valued superfluid ground states^’^. These ap¬ 
proaches either discard correlations'^ or are restricted to flnite- 
size systems and the standard Bose-Hubbard model^. In a 
recent work on a double well using a mean-fleld model a 
complex ground state has been found for large pair-tunneling 
amplitudes'^. Our theoretical framework allows a correlated 
treatment of a lattice in the thermodynamic limit incorporat¬ 
ing beyond-Hubbard processes and offers a conclusive under¬ 
standing of the twisted superfluid quantum phase. 

At first glance, the complex order-parameter and the time- 
reversal symmetry breaking seem to stand in contradiction to 
Feynman’s no-node theorem^’^, which we resolve in the fol¬ 
lowing. For real Hamiltonians, the no-node theorem implies a 
positive real and therefore non-degenerate ground state. How¬ 
ever, the energy difference between the lowest-energetic states 
can become arbitrarily small with increasing system size. It 
has been shown that the energy gap between the lowest two 
states vanishes in the thermodynamic limit for the charge den¬ 
sity wave^^ as well as for a highly occupied double well^^. In 
the first case, the superposition of two opposite density waves 
represents the ground state in finite systems but a small per¬ 
turbation breaks the inversion symmetry in favor of one den¬ 
sity wave. In the second case, the time-reversal symmetry 
becomes unstable against perturbations for a large number of 
particles. The mean-fleld approach reproduces the thermody¬ 
namic limit and tends to break both symmetries via the inco¬ 
herent mean-fleld coupling at the boundary of the cluster. For 
real macroscopic systems, such incoherences are introduced 
by several mechanisms such as system-bath coupling with the 
surrounding, thermodynamic excitations, fragmentation, infi¬ 
nite relaxation times and dissipation in open systems. 

The driving mechanism of the twisted superfluid phase is 
the correlated pair tunneling P competing with the single¬ 
particle tunneling J. The competition becomes apparent when 
introducing simple variational mean-fleld expressions. Re¬ 
placing all operators with complex numbers (see 

Methods), we can approximate both contributions as 

— J(dtdj) + c.c. ~ — 2J(j)i(l)j cos(d), 

P/2 (dtdtdjdj) + c.c. ^ P(l)^(t>^j cos(2d) 

with 0 = 0i—0j. For small amplitudes of the pair tunneling 
P, the first term dominates and minimizing the energy yields 
d = 0. The pair tunneling energy on the other hand is mini¬ 
mized for 0 = ±7r/2. Due to the different functional depen¬ 
dence on 0, the system undergoes a phase transition from 
0 = 0 to a finite value |6>|>0 at a critical ratio P/J. In a 
more qualitative picture, we have to include correlations that 
are strong near the insulating phases, i.e. {hi) and {hihi) may 
differ strongly from the variational mean-fleld expressions 
and respectively. This is accounted for in the applied 

cluster Gutzwiller approach treating the cluster sites in a full 
many-particle description. Additionally, the off-site interac¬ 
tion V creates a competition between the phase twist and den- 
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FIG. 2: Phase boundaries and order parameters, a Cut through 
the phase diagram Fig. lb showing the plane with P = V = 0.15U 
(see Fig. lb for abbreviations). The markers indicate the critical point 
for the SS-TSS phase transition obtained for cluster sizes 1 (triangle), 
2 (circle) and 4 (square), b-d Details of the phase-diagram showing 
dimerized density-wave phases (DDW) with fillings 1/4, 3/4 and 
5/4 as well as supersolids with long periodicity (DSS). e Averaged 
superfluid order parameter | {a) \ (black), its complex phase angle 0 
(yellow), density wave | An| (blue), and the energy derivative (red) 
shown across various phase transitions for J/U = 0 . 12 , where 7 is 
ground-state degeneracy. 


sity modulations. The interplay of insulators, density waves 
and complex twisting of the order parameter is vital to the 
phase diagram in Fig. lb. Asa result, a multitude of quantum 
phase transitions appear which we discuss in the following. 

Figure 2a shows a cut through the phase diagram for 
P = V = O.lbU, which depicts also the supersolid phase (SS) 
omitted in Fig. lb. The supersolid has a non-vanishing su¬ 
perfluid order parameter combined with a density wave rep¬ 
resenting diagonal order^^’^^“^^’^^“'^^(Fig. Ig). The respective 
expectation values for (d^) and An = {hi — hj) are shown for 
J = 0.12 in Fig. 2e. In the thermodynamic limit, the super¬ 
solid phase breaks the inversion symmetry P and has a two¬ 
fold degenerate ground state { 7 /, PV^}, where the operation P 
interchanges the two hexagonal sublattices A and B. 

The transition between the supersolid and the twisted super¬ 
fluid occurs via a twisted supersolid (TSS) which combines 
the properties of both phases. In the twisted supersolid both 
the density wave and the complex twisting d 7 ^ 0 of the order 
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parameter occur (Fig. If). Both inversion and time-reversal 
symmetry are broken and the ground state is four-fold de¬ 
generate { 7 /;, Pt/;, Tt/;, PTt/)}. The superfluid order parame¬ 
ter, density wave and phase twist 0 have a kink at the SS-TSS 
as well as at the TSS-TSF transition (Fig. 2e). In the TSS 
phase the phase twist 0 increases rapidly towards ±7r/2 with 
increasing /ieff, whereas the density modulation drops to zero 
when entering the TSF phase. This means that the complex 
phase of the order parameter suppresses the density wave re¬ 
flecting the interplay between nearest-neighbor interaction V 
and pair tunneling P. The energy derivative —dH/djj. has 
no discontinuities but kinks at transitions between SF and SS 
phases marking all phase transitions as second (or higher) or¬ 
der. For the transition from the TSS to the SS, we can evaluate 
the energy gap A between the ground state and the real-valued 
SS state. At the phase boundary, we And the critical behavior 
A (X l/i - \ J - and |P - Pd^^ with a critical 

exponent zu = 1.96 ± 0.10. We observe that the TSS phase 
boundary expands when increasing the cluster size of the cal¬ 
culation (Fig. 2a) indicating that flnite-size scaling"^^ would 
lead to an even larger region where the superfluid order param¬ 
eter is complex valued. Note that Fig. lb is a zero-temperature 
phase diagram. We expect the energetically lowest excitation 
within the twisted superfluid phase to be a disturbance of the 
phase angle, which would be associated with an energy on the 
order of the pair tunneling P. 

Dimerized density-wave insulators 

For small J/U several insulator phases can be found with a 
vanishing order parameter (d) (Figs. 1 and 2). For V = P = 0, 
we And the usual Mott insulators (MI)^^“^^ with integer Ailing 
(nd = 1, 2, 3 on A and B sites. For increasing off-site interac¬ 
tion, charge-density-wave insulators (CDW)^^’^^“^^ with half¬ 
integer Ailing per lattice site appear in between the Mott lobes. 
These incompressible phases are characterized by an alternat¬ 
ing Ailing n and n + 1 on A and B sites (Fig. Ic). Like in 
the supersolid phases, the inversion symmetry of the lattice is 
spontaneously broken in the thermodynamic limit. Respon¬ 
sible for this is the nearest-neighbor interaction V which is 
reduced for an alternating Ailing. The charge-density wave 
phases are also referred to as Mott solids, alternating Mott in¬ 
sulators, and for square lattices as checkerboard insulators. 

Remarkably, in between charge-density waves and Mott in¬ 
sulators, we also And fractional insulator phases with quarter- 
integer Ailing indicated in blue in Figs. 1 and 2. These 
dimerized density waves (DDW) combine a density wave 
and a vanishing superfluid order parameter with a spatial 
delocalization of particles on dimers of the lattice. Simi¬ 
lar phases have been discussed only for lattices with non- 
symmetric tunneling and are referred to as ’’loophole” or frac¬ 
tional insulators'^^"'^^ or in the context of honeycomb lattices 
as dimerized Mott insulators^^. For translational-symmetric 
lattices this kind of dimerization without a density wave is en¬ 
tirely unexpected. We And that for non-zero nearest-neighbor 
interaction V dimerized phases with quarter-integer Ailing can 
emerge from the spontaneous breaking of the translational 


symmetry. For a Ailing of 1/4, a long-range density wave 
emerges with a delocalized particle on two sites alternating 
with two empty sites (illustrated in Fig. li). On these four 
sites, the respective many-particle state can be written approx¬ 
imately as (|0,1,0,0) + |1,0,0,0) )/\/2. Due to six possi¬ 
ble conflgurations on the honeycomb lattice, this state is six¬ 
fold degenerate. Like the CDW, the dimerized density-wave 
minimizes the V energy but also allows the particles gaining 
energy from the single-particle tunneling J. The flrst-order 
energy L^ddWi/4 = —J — per unit cell can fall below the 
energy of the CDW L^cdWio = — 2/i. 

The close-ups of the DDW phases for Ailing 1/4, 3/4 and 
5/4 in Fig. 2b-d reveal the typical droplet shape that is also 
present for non-symmetric lattices^^"'^^ as well as a direct 
insulator-insulator transition in Fig. 2b. We also And this kind 
of dimerization in the supersolid state. Close to the dimerized 
density-wave with Ailing 1/4 and 3/4, a phase emerges with 
a non-vanishing superfluid order parameter and equal densi¬ 
ties on neighboring dimer sites but a long-range density wave 
(Fig. 2b,c). It is denoted as dimerized supersolid (DSS) and 
illustrated in Fig. Ih. 

In the dimerized density-wave, neighboring sites are 
strongly correlated and thus single-site mean-fleld approaches 
are unable to And this quantum phase. Due to the cluster size 
for this calculation. Ailing at multiples of n = 1/4 are ex¬ 
pected, but depending on the lattice geometry, different flll- 
ings might be possible and energetically favorable depending 
on the chemical potential. For example, the structure shown in 
Fig. Ij has the same energy per particle as the DDW 1/4 state 
but with a larger Ailing 5/18 indicating that fragmentation into 
larger unit cells might be possible. 

Before turning to two-component systems, we briefly dis¬ 
cuss other phases appearing in Figs. 1 and 2 that are already 
discussed in the literature. At a critical value V = U/6 it 
becomes energetically favorable to gather all particles on one 
sublattice, while for smaller values of V the population imbal¬ 
ance between the sites is minimized^"^. As a result, both Mott 
and CDW insulators undergo the transition to a maximally- 
imbalanced charge-density wave (CDW2). For large values 
of V, the dimerized density waves (with Ailing >3/4) have a 
smooth crossover to a long-range density-wave with maximal 
imbalance between neighboring sites, i.e. |n, 0, n — 1,0). For 
sufflciently large values of //eff the insulators develop a flnite 
pair superfluidity (d^) driven by the pair tunneling P. Fig¬ 
ure lb shows that pair superfluids (PSF)^^“^^ extend from the 
Mott insulator phases, whereas CDW insulators are attached 
to pair supersolids characterized by both a den¬ 

sity wave and pair superfluidity. 

Twisted superfluidity for two components 

In the following, we discuss the quantum phases described 
above for a system with two distinguishable bosonic species. 
This introduces a new degree of complexity and has a dra¬ 
matic impact on the phase diagram. The two bosonic compo¬ 
nents, denoted as t and can e.g. be realized as atoms in two 
different hyperflne states. For simplicity, we assume symmet- 
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FIG. 3: Phase diagrams of the extended two-component Bose- 
Hubbard model. The results are obtained for a P = 0.02 U and b 
P = 0.05 U with P = V = C using a single-site cluster. The close- 
ups below for c P = 0.02 U and d P = 0.05 P are computed with 
a two-site cluster. This allows identifying dimerized density waves 
(DDW) and dimerized supersolids (DSS). e Order parameters and 
degeneracy 7 for J = 0.1 P and P = 0.05 P. In the continuously 
degenerate phases the expectation values correspond to the solution 

with 0 = 6>max. 


ric intra- and interspecies interactions as well as an inversion- 
symmetric honeycomb lattice. Note that in the experiment^ 
the honeycomb lattice is strongly spin-dependent with a site- 
offset between the sublattices. Furthermore, we assume a bal¬ 
anced population of both species, i.e 
The two-component Hamiltonian is given by 

H = H^ + Hi + uJ2 ’ (3) 

i 

where and are the intraspecies Hamiltonians (1). The 
nearest-neighbor interspecies interaction reads 

H ^ a\ialiaijafj+C (4) 

{'id) 

and accounts for the nearest-neighbor interaction V, the tun¬ 
neling of a t, i-pair with amplitude P and the cross tunnel¬ 
ing (or counter hopping) C, which interchanges a t- and a 
l-particle on neighboring sites. For spin-independent contact 
interaction, the amplitudes V = P = C are equal to the cor¬ 
responding intraspecies interactions. The complexity of the 
two-component Hamiltonian is reflected by eight independent 


mean-fleld parameters that are needed for the description (see 
Methods). In addition, the number of the many-particle states 
in the cluster increases drastically by incorporating the tensor 
basis of the two components limiting the cluster size to maxi¬ 
mally two sites. 

The phase diagram for the binary bosonic mixture is shown 
in Fig. 3 and differs strongly from the single-component case. 
Surprisingly, the twisted superfluid phase (TSF) appears al¬ 
ready for small values of the chemical potential and, most 
notably, for small amplitudes of the pair tunneling. For 
P = 0.05(7, the TSF phases form stripes in the phase diagram 
emerging from the charge density waves (CDW) (Fig. 3b) and 
are separated by supersolid phases (SS). For a smaller am¬ 
plitude P = 0.02(7 the stripes overlap as shown in Fig. 3a. 
We observe that the slope of the SF-TSF phase boundary is 
roughly proportional to P. To first order, this correspond to a 
constant ratio (P)/(J) ^ P{h)‘^/J{h) leading to J oc (n)P. 

The twisted superfluid phase for two components is char¬ 
acterized by anti-aligned twisted complex order parameters 
with 0^ + = 0 as illustrated in Fig. 4a,b. The state is 

spontaneously symmetry-broken and forms degenerate pairs 
with 0^ = —0^ = P0. This implies that the two components 
are in a strongly correlated state with respect to each other. 
The phase correlation is established by the cross tunneling 
process (equation (4)). Using the simple mean-fleld ap¬ 
proach (2), the cross-tunneling energy is approximated by 
2C(l)^i(t)lj(t)ii(t)^j cos {0^ — which becomes minimal for 
0^ = —0l = 7r/2. The actual values of 0^ and are de¬ 
termined by all contributions in (3) and (4). The cross tun¬ 
neling energy for a correlated state C(d|-d|^ d^id^j) can be 
large even in the vicinity of the Mott phases. In the pure 
mean-fleld approach of Ref. 4, it is pointed out that this 
beyond-mean-fleld effect could possibly yield a twisted su¬ 
perfluid phase. Our results show that the cross tunneling pro¬ 
cess is indeed mainly responsible for the emergence of the 
TSF phase at small amplitudes P. At the same time so-called 
spin-density waves may appear, which represent anti-aligned 
single-component density waves with An^ = —An^ such that 
the total density is homogeneous An^pAni = 0. 

The experimental observation of the density wave requires 
imaging with single-site resolution, whereas the complex 
twisting of the order parameter can be revealed in momentum 
space. Time-of-flight experiments with ultracold atoms allow 
the mapping on free momenta after a sudden release of the 
atomic trapping potential. For the individual spin components 
(we omit the spin index), the momentum density is propor¬ 
tional to the structure factor"^^ S'(k) = (dtdj), 

which gives six equivalent first-order momentum peaks Sq 
at reciprocal lattice vectors for a superfluid (0 = 0). In 
the twisted superfluid, this symmetry is broken and adjacent 
peaks have either amplitude 5'+ or S-. Figures 4c,d de¬ 
pict calculated time-of-flight images yielding alternating pat¬ 
terns in momentum space. The two amplitudes are given by 
S±{0(j) = 25'o[l + ^cr cos(^cr ± 7r/3)] for a either f or 
For two components, the orientation of the stronger peaks are 
opposite for both components due to S±{0^) = Szf{0i). This 
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FIG. 4: Experimental signature of twisted superfluidity. For two 

components, denoted as t and the twisted superfluid phase is char¬ 
acterized by an anti-correlation of phase twist and density wave, 
a The arrows illustrate a complex phase angle 0^ = 0 for the t- 
component. b The ^-component has the opposite twist = —0. The 
density waves oppose each other forming a spin-density wave with 
a homogeneous total density. c,d Illustration of time-of-flight im¬ 
ages allowing the experimental observation of the twisted superfluid 
phase. The images correspond to the individual spin components of 
a and b, respectively. The alternating intensities of the first-order 
coherence peaks interchange between the species due to 0^ = —0^. 
Here we chose O^tt/S and = 1, whereas at ^ = tt/S the weaker 
peaks disappear entirely, e The relation between phase twist 0 and 
spin-density wave for the infinitely degenerate ground state in the 
TSF phase. The markers correspond to 200 individual solutions of 
the self-consistent cluster mean-field algorithm with different ran¬ 
dom initial states. 


behavior was also observed experimentally^. An additional 
spin-density wave An^ modifies the strengths of the peaks 
via with whereas A^j = 1 for a homogeneous density (see 
Methods). 

Infinite degeneracy and quantum phase transitions 

In contrast to the two-fold degeneracy for a single species, 
we find that the TSF ground state for two components has an 
infinite degeneracy in respect to the continuous parameter 0. 
More precisely, all states with — ^max ^ 0 < 6>max have the 
same energy, where 0 can take a continuous value and ^max 
depends on the point in the phase diagram. As depicted in 
Fig. 4e, the spin-density wave An^ = — An^ and the complex 
twisting 0 are correlated. This corresponds to a redistribu¬ 
tion of the energy between the density wave and the complex 
twisting of the superfluid order parameter. For 0 = ±^max the 
spin-density wave vanishes, whereas the spin-density wave 
I Aricrl is maximal for 6> = 0. These two extrema correspond 
to a ’’pure” twisted superfluid with Aua- = 0 and a real, anti¬ 
aligned supersolid in each component for 6> = 0, respectively. 
The latter implies that two of the solutions are real wave func¬ 
tions obeying time-reversal symmetry. The origin of the con¬ 
tinuous symmetry lies within an equal redistribution of the en¬ 
ergy between the different interaction processes and persists 
only if all amplitudes V, P and C are equal, which holds for 
equivalent atomic species with SU(2) symmetric interactions. 
We observe this degeneracy for all studied cluster sizes. How¬ 
ever, it might be lifted by including higher bands of the lattice 


leading to a renormalization of the amplitudes above^^. At the 
touching points of TSF phases and CDW insulators, we ob¬ 
serve a complex order parameter in combination with a total 
density wave An^ + An^^ ^ 0, i.e. a twisted supersolid (TSS). 
In analogy to the TSF phase, its ground state is continuously 
symmetry broken. We cannot exclude the possibility of a thin 
twisted supersolid layer at the boundary between supersolid 
and TSF phase. 

Similar, we find also an infinite degeneracy of the ground 
state for the MI and CDW insulator phases for symmetric 
intra- and interspecies interaction. For instance, the Mott in¬ 
sulator with filling {n^i -h = 1 on two neighboring sites 
i and j can be written in first order as (ci |t) • + C 2 |i) J C) 
(c2 |t)j ~ Cl 11)^ ) with an arbitrary constant ci and |cip + 
|c 2 p = 1. Here, the intra- and interspecies interaction en¬ 
ergy of nearest-neighbor sites V and the cross tunneling en¬ 
ergy C are redistributed yielding the same total energy. Since 
the cross-tunneling order parameter (d|-d^i) does not vanish, 
the two-species Mott insulator is sometimes referred to as a 
counterflow superfluid 

As shown in Fig. 3e, the transitions between SS and TSF 
phases are first order quantum phase transitions with discon¬ 
tinuities in d{H)/dji. At the transition point the total density 
wave An^ + An^ and the superfluid order parameter are dis¬ 
continuous. The structure of the ground state changes across 
the transition from a two-fold to infinite degeneracy. Also the 
transition between superfluid and Mott insulator are of first 
order^^’^^ characterized by a discontinuity of the superfluid or¬ 
der parameter as well as by the different degeneracies. 

The close-up of the phase diagram shown in Fig. 3d reveals 
the appearance of dimerized density waves (DDW) for two 
components. These phases with quarter-integer filling arise 
between Mott and CDW insulators as discussed above for a 
single component. Due to the strong correlations between 
neighboring sites, this quantum phase is only found in cal¬ 
culations using clusters with at least two sites. Similar as the 
Mott insulator phase, the DDW has an infinite degeneracy. 
Attached to this phase, we also find a dimerized supersolid 
(DSS) in analogy to the single-component phase diagram. 

Discussion 

The appearance of the twisted superfluid phase at small val¬ 
ues of the pair tunneling amplitude for a bosonic mixture 
is quite remarkable. For the honeycomb lattice a pair tun¬ 
neling amplitude of P/U = 0.05 and 0.01 correspond to a 
lattice depth Vq = to 3.5 F^r, respectively^"^. The ra¬ 
tio J/U depends on the scattering length a^, where we find 
J/U = 0.14 for Vo = 3.5F^r and = 360ao. However, 
the ratios P/J and C/J can also be tuned experimentally us¬ 
ing lattice shaking techniques'^. The presented results shed 
light on the formation of twisted superfluidity recently found 
experimentally^. We are able to draw a consistent picture of 
the time-reversal symmetry breaking and the interplay with 
the insulating phases. The strikingly complex phase diagram 
is a manifestation of the multilayered interplay between differ¬ 
ent extended Hubbard processes on the one hand and nearest- 
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neighbor correlations on the other hand. In the experimental 
realization^, where the scattering length is about = lOOao, 
the lattice potential depends on the atomic hyperfine states 
causing an energy offset between the sublattices. In future 
studies, the effect of the spin-dependent offset on the phase 
diagram can be evaluated. 

Spontaneous breaking of the translational symmetry allows 
for dimerized density waves with quarter-integer filling. The 
nature of this fractional insulator is reflected in a strongly cor¬ 
related superposition state on dimers of the lattice, which is 
linked to resonating valance bond coupling in solids. The 
same mechanism leads to dimerized supersolid phases. We 
have not found any indication for phase separation in the hon¬ 
eycomb lattice similar as Refs. 30,31. 

It would be interesting to study in detail the infinite degen¬ 
eracy for the two-component TSF and TSS ground states. This 
degeneracy should be experimentally observable for a spin- 
independent lattice and SU{N) symmetric interactions. For 
the hyperfine states of rubidiums atoms^, the SU (N) symme¬ 
try is broken causing slightly different intra- and interspecies 
scattering length. This also breaks the continuous ground- 
state degeneracy of the TSF phase and locks the value of 0 
and the density wave. However, the two-fold degeneracy con¬ 
nected with the time-reversal symmetry breaking persists. 

For the sake of simplicity and to avoid an even larger pa¬ 
rameter space, the Hamiltonians (1) and (3) are not complete 
with respect to nearest-neighbor interaction processes. They 
neglect a term of the form a]{hi + hj)ai known as density- 
induced tunneling. For ultracold neutral atoms, its ampli¬ 
tude is comparably large. However, it has been shown that 
it may be absorbed in the single-particle tunneling J in most 
situationsThus, we mainly expect a shift of the phase 
boundaries due to this combined occupation-dependent tun¬ 
neling. Moreover, the next-nearest-neighbor tunneling has a 
considerable amplitude in shallow lattices^'^. Since its sign 
corresponds to antiferromagnetic coupling, it may even sta¬ 
bilize the twisted superfiuid phase. Long-range interactions, 
such as dipolar interaction^^, offer further perspectives for 
studying off-site interactions. 

The twisted superfiuid itself presents a highly exotic state 
of matter as it combines the off-diagonal long-range order of 
a superfiuid state with a correlated short-range ordering. For 
two components, correlations between the species are essen¬ 
tial and manifest themselves in spin-density waves and an op¬ 
posite complex twisting of the superfiuid order parameter. As 
previously mentioned the connection to time-reversal symme¬ 
try breaking in the context of superconductors^’is of high 
interest. Quantum gas experiments may contribute by offering 
highly tunable systems and a different experimental platform 
for studying correlated quantum phases in lattices. 

Methods 

Cluster mean-field method and order parameters. For the com¬ 
putation of the phase diagram, we apply the cluster Gutzwiller 
method. It is based on factorizing the many-body wave func¬ 
tion 1^) l^c) in two parts describing particles inside |^c) and out¬ 
side the cluster |^). Accordingly, the Hamiltonian decomposes in 


Hnm = {M\ iFciuster + (('0| ^boundary |V^)) |iV) in thc basis of many- 
particle cluster states \N). The knowledge of the (factorized) solu¬ 
tion for 1 ^) allows obtaining the correlated cluster wave function 
l^c) = EiV cn\N) by means of diagonalization. Fortunately, 
tight-binding Hamiltonians reduce the information required from the 
wave function |^) outside the cluster, since ^boundary acts only on 
sites at the boundary. For the single-component Hamiltonian (1), 
^boundary \'^) rcduccs to (tti), {hi) and {didi), which we obtain 
self-consistently from inside the cluster. Note that ^boundary couples 
different Fock spaces in the cluster. In the case of the two-component 
Hamiltonian (3) eight mean-field parameter are required, i.e. (dfi), 
(ntz), (d^id^i), (dii), (%), (diidii), (d^id^i) and {d\.dii). We 
apply self-consistent iterations /, until \cn,i — cn,i-2 \ falls be¬ 
low a threshold value (typically smaller than 10“®). On subsequent 
iterations, we allow the cluster wave function and mean fields to be 
different, which enables us to find structures spanning 2k sites with a 
/c-site cluster. Due to the high dimension of the mean-field parameter 
space up to several thousand iterations are needed in certain regions 
of the phase diagram. We have carefully checked convergence with 
respect to the size of the cluster basis. Further technical details can 
be found in Ref. 41. 

Phase diagrams. The phase diagrams presented in Figs. 1 and 2 
are obtained by means of the cluster mean-field approach using a 
two-site cluster. The phase boundaries are slightly affected when 
increasing the cluster size (see comparison in Fig. 2a), since longer- 
ranged many-particle correlations can be incorporated^ k The effec¬ 
tive chemical potential //eff = (/i + f//2)/(l -h 6P/U) — U/2 ac¬ 
counts for the increased total interaction energy for finite V. In these 
units, the position of the Mott lobes is approximately independent of 
P = V. The phase diagram for two components in Fig. 3a,b is com¬ 
puted using a single-site cluster and the two-component tensor basis, 
i.e. including all interspecies on-site correlations. The close-ups in 
Fig. 3c,d are obtained using two-sites clusters. For two components, 
the chemical potential /Xeff = {/j. U/2)/ {1 3P/U) — f//2 ac¬ 
counts for the off-site interaction energy. 

Simple variational mean-field expressions. For a single compo¬ 
nent, a variational model can be constructed serving as a motiva¬ 
tion for the complex phase in the twisted superfiuid phase. In con¬ 
trast to the cluster mean-field approach above, we neglect all cor¬ 
relations arising from the one- and two-particle operators by set¬ 
ting (dldj) = {di)*(dj) and (dld^jdkdi) = (di)*(dj)*{dk)(di) . Al¬ 
lowing for two different sublattices A and B, we replace the oper¬ 
ators with complex numbers (aA) = 0 a and (ds) = 0 Be^^ with a 
relative phase 0. In this model, the mean-field energy of Hamilto¬ 
nian (1) per unit cell reads {H) = L(0a+0b)/2 — m(0a+0b) “ 
6J0A0B cos(^) + 6y0A0B + 3H0A0B cos(2^). This model pre¬ 
dicts the transition between superfiuid and twisted superfiuid qual¬ 
itatively, but is unable to predict correlated insulator phases and is 
invalid for small J/U. The same model can also be applied to the 
two-species case (equation (3)) with the same restrictions mentioned 
above. In particular, this approach cannot capture the twisted super¬ 
fiuid phase at small amplitudes of the off-site interactions, which is 
driven by correlations. 

Structure factor. The twisted superfiuid state gives rise to a distinct 
time-of-fiight signature showing first-order momentum peaks with 
an alternating intensity S± (Fig. 4). A density wave Aua within 
the spin component a modifies the strengths of the peaks S± via 
Act = (1 — with N = + haj). The central 

peak for zero momentum is proportional to 1 + Act cos(^i). 
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